import numpy as np
import os
# local modules
import geom
geom = reload( geom )
import abaqus_functions
abaqus_functions = reload( abaqus_functions )
# abaqus functions
from abaqus import *
from regionToolset import Region
# local constants
from constants import *
# abaqus constants
from abaqusConstants import *
#
def create_model( self ):
    # abaqus modules
    self.rebuild()
    R          = self.r # mm
    H          = self.h # mm
    R2         = self.r2 # mm
    PLYT       = self.plyt # mm
    LAMINAPROP = self.laminaprop
    STACK      = self.stack
    R_NUMEL    = self.numel_r
    ELSIZE_R   = self.elsize_r # mm
    ELSIZE_H   = self.elsize_h # mm
    AXIALDISPL = self.axialdispl # mm
    JOBNAME    = self.jobname
    #
    backwardCompatibility.setValues( includeDeprecated = True,
                                     reportDeprecated  = False )
    # 
    #Mdb()
    #mdb.saveAs(pathName=r'C:\Temp\abaqus\cylinder.cae')
    self.mdb = mdb
    #RESETING THE CURRENT VIEWPORT
    myView = session.viewports[ session.currentViewportName ]
    myView.setValues( displayedObject = None )
    #CREATING A NEW MODEL
    mod = mdb.Model(name=JOBNAME)
    self.mod = mod
    #DELETING THE DEFAULT MODEL
    tmp = 'Model-1'
    if tmp in mdb.models.keys():
        del mdb.models[tmp]
    #CREATING A NEW PART
    part = mod.Part( name = 'Cylinder',
                     dimensionality = THREE_D, 
                     type = DEFORMABLE_BODY )
    self.part = part
    #CREATING THE SKETCH which will be used to create the shell geometry
    s1 = mod.ConstrainedSketch( name='SketchCylinder', 
                                sheetSize=max( [2.1*R, 1.1*H] ) )
    #sketch profile to be extruded
    s1.CircleByCenterPerimeter( center = (0.0, 0.0),
                                point1 = (  R, 0.0) )
    #
    #CREATING A LOCAL COORDINATE SYSTEM TO USE IN THE BOUNDARY CONDITIONS
    csysCyl = part.DatumCsysByThreePoints( name='CSYSCylinder',
                                           coordSysType=CYLINDRICAL,
                                           origin=(0,0,0),
                                           point1=(1,0,0),
                                           point2=(1,1,0) )
    #CREATING THE CYLINDER SHELL GEOMETRY
    myCyl = part.BaseShellExtrude( 
                                   depth = H, 
                                   draftAngle = -self.alphadeg,
                                   sketch=s1,
                                 )
    #CREATING PARTITIONS considering holes and ploads
    thetas = []
    for hole in self.holes:
        thetas.append( hole.theta1 )
        thetas.append( hole.theta2 )
    for pload in self.ploadconf.ploads:
        thetas.append( pload.theta )
    # the planes will cut both sides, so the thetas higher than 180.
    # must be cut at the opposite side
    thetas2 = set() 
    for theta in thetas:
        if theta >= 180.:
            thetas2.add( theta - 180. )
        else:
            thetas2.add( theta )
    thetas2 = list(thetas2)
    thetas2.sort()
    # creating points and instantiating cutting plane objects
    part.DatumPointByCoordinate( coords = (0.0, 0.0,   0.0) ) 
    po = part.datums[3]
    part.DatumPointByCoordinate( coords = (0.0, 0.0, 1.1*H) )
    pz = part.datums[4]
    i = 0
    thetas = []
    # adding +180. to all thetas in order to create a theta vector
    # which will allow correct edge and face identification
    self.cutplanes = []
    for theta in thetas2:
        i += 1
        x = 1.1*R * np.cos( np.deg2rad( theta ) )
        y = 1.1*R * np.sin( np.deg2rad( theta ) )
        part.DatumPointByCoordinate( coords = (x, y, 0.0) )
        pp = part.datums[ 4 + i ]
        cutplane = geom.CutPlane()
        cutplane.p1 = po
        cutplane.p2 = pz
        cutplane.p3 = pp
        cutplane.part = part
        self.cutplanes.append( cutplane )
        thetas.append( theta      )
        thetas.append( theta+180. )
    # creating cutting planes
    for cutplane in self.cutplanes:
        cutplane.create()
    # cutting faces     
    for cutplane in self.cutplanes:
        tmp = part.PartitionFaceByDatumPlane( datumPlane=cutplane.datum,
                                              faces=part.faces )
    # dividing edges for ploads
    rmid, zmid   = self.r_z_from_ppoint( 0.05 )
    rmid2, zmid2 = self.r_z_from_ppoint( 0.95 )
    rbot, zbot   = self.r_z_from_ppoint( 0.00 )
    rtop, ztop   = self.r_z_from_ppoint( 1.00 )
    #
    ploadsthetas = []
    for pload in self.ploadconf.ploads:
        theta = pload.theta
        ploadsthetas.append( theta )
        x  = rmid * np.cos( np.deg2rad( theta ) )
        y  = rmid * np.sin( np.deg2rad( theta ) )
        z  = zmid
        edge = part.edges.findAt( ((x, y, z),) )
        if pload.theta < 180.:
            part.PartitionEdgeByParam( edges = edge,
                                       parameter = ( 1. - pload.pt) ) 
        else:
            part.PartitionEdgeByParam( edges = edge,
                                       parameter = ( pload.pt     ) ) 

    # part edges identification
    self.part_edges = []
    for theta in thetas:
        edge = geom.Edge()
        edge.theta = theta
        x = rmid * np.cos( np.deg2rad( theta ) )
        y = rmid * np.sin( np.deg2rad( theta ) )
        z = zmid
        edge.mid = part.edges.findAt( ((x, y, z),) )
        if theta in ploadsthetas:
            x = rmid2 * np.cos( np.deg2rad( theta ) )
            y = rmid2 * np.sin( np.deg2rad( theta ) )
            z = zmid2
            edge.mid2 = part.edges.findAt( ((x, y, z),) )
        x = rbot * np.cos( np.deg2rad( theta + 2. ) )
        y = rbot * np.sin( np.deg2rad( theta + 2. ) )
        z = zbot
        edge.bot =  part.edges.findAt( ((x, y, z),) )
        x = rtop * np.cos( np.deg2rad( theta + 2. ) )
        y = rtop * np.sin( np.deg2rad( theta + 2. ) )
        z = ztop
        edge.top = part.edges.findAt( ((x, y, z),) )
        #
        self.part_edges.append( edge )
    self.part_edges.sort()
    #FACES identification
    self.faces = []
    for theta in thetas:
        x = rmid * np.cos( np.deg2rad( theta + 2. ) )
        y = rmid * np.sin( np.deg2rad( theta + 2. ) )
        z = zmid
        face = geom.Face()
        face.theta = theta
        face.feature = part.faces.findAt( ((x, y, z),) )
        self.faces.append( face )
    self.faces.sort()
    #_________________________________________________________
    #
    # lamina
    #_________________________________________________________
    #
    #CREATING A LAMINA MATERIAL
    localCsys = part.datums[1]
    if self.using_lam_parameters == False:
        myMat = mod.Material( name='carbonfiber' )
        #                      
        myMat.Elastic( table = ( LAMINAPROP , ) ,
                       type  = LAMINA )
        ALLOWABLES = tuple([ abs(i) for i in self.allowables ]) 
        myMat.HashinDamageInitiation( table = ( ALLOWABLES, ) )
        #CREATING A COMPOSITE LAYUP
        myLayup = part.CompositeLayup(\
                          name='CompositePlate',
                          description='stack from inside to outside',
                          offsetType = MIDDLE_SURFACE,
                          symmetric = False,
                          thicknessAssignment = FROM_SECTION,
                          elementType = SHELL,
                        )
        myLayup.Section(\
                         preIntegrate = OFF,
                         integrationRule = SIMPSON,
                         thicknessType = UNIFORM,
                         poissonDefinition = DEFAULT,
                         temperature = GRADIENT,
                         useDensity = OFF,
                       )
        myLayup.ReferenceOrientation(\
                         orientationType = SYSTEM,
                         localCsys = localCsys,
                         fieldName = '',
                         additionalRotationType = ROTATION_NONE,
                         angle = 0.,
                         additionalRotationField = '',
                         axis = AXIS_2,
                       )
        #CREATING ALL PLIES
        icount = 0
        for face in self.faces:
            icount += 1
            jcount = 0
            rFace = Region( faces = face.feature )
            for angle in STACK:
                jcount += 1
                myLayup.CompositePly(\
                     suppressed = False,
                     plyName = 'face_' + str(icount) +\
                               '_ply_' + str(jcount),
                     region = rFace,
                     material = 'carbonfiber',
                     thicknessType = SPECIFY_THICKNESS,
                     thickness = PLYT,
                     orientationValue = angle,
                     orientationType = SPECIFY_ORIENT,
                     numIntPoints = 3,
                    )
    #DEFINING THE MESH SEEDS ALONG ALL EDGES
    for edge in self.part_edges:
        part.seedEdgeBySize( edges = edge.mid,
                             size = ELSIZE_H,
                             constraint = FREE )
        part.seedEdgeBySize( edges = edge.bot,
                             size = ELSIZE_R,
                             constraint = FREE )
        part.seedEdgeBySize( edges = edge.top,
                             size = ELSIZE_R,
                             constraint = FREE )
        if edge.mid2 <> None:
            part.seedEdgeBySize( edges = edge.mid2,
                                 size = ELSIZE_H,
                                 constraint = FREE )
        
    # defining mesh controls
    part.setMeshControls( regions = part.faces,
                          elemShape = QUAD,
                          technique = STRUCTURED )
    import mesh
    if self.elem_type == 'S4R':
        elemtype = mesh.ElemType( elemCode = S4R,
                                  elemLibrary = STANDARD )
    elif self.elem_type == 'S8R':
        elemtype = mesh.ElemType( elemCode = S8R,
                                  elemLibrary = STANDARD )
    elif self.elem_type == 'S8R5':
        elemtype = mesh.ElemType( elemCode = S8R5,
                                  elemLibrary = STANDARD )
    else:
        elemtype = mesh.ElemType( elemCode = S4R5,
                                  elemLibrary = STANDARD )
    part.setElementType( regions = Region(faces=part.faces),
                         elemTypes = (elemtype,) )
    # assemblying the cylinder - creating instance
    inst = mod.rootAssembly.Instance( name='InstanceCylinder',
                                      part=part,
                                      dependent=ON )
    # finding instance edges
    self.inst_edges = []
    for theta in thetas:
        edge = geom.Edge()
        edge.theta = theta
        x = rmid * np.cos( np.deg2rad( theta ) )
        y = rmid * np.sin( np.deg2rad( theta ) )
        z = zmid
        edge.mid = inst.edges.findAt( ((x, y, z),) )
        if theta in ploadsthetas:
            x = rmid2 * np.cos( np.deg2rad( theta ) )
            y = rmid2 * np.sin( np.deg2rad( theta ) )
            z = zmid2
            edge.mid2 = inst.edges.findAt( ((x, y, z),) )
        x = rbot * np.cos( np.deg2rad( theta + 2. ) )
        y = rbot * np.sin( np.deg2rad( theta + 2. ) )
        z = zbot
        edge.bot =  inst.edges.findAt( ((x, y, z),) )
        x = rtop * np.cos( np.deg2rad( theta + 2. ) )
        y = rtop * np.sin( np.deg2rad( theta + 2. ) )
        z = ztop
        edge.top = inst.edges.findAt( ((x, y, z),) )
        #
        self.inst_edges.append( edge )
    self.inst_edges.sort()
    # finding vertices for ploads
    for pload in self.ploadconf.ploads:
        pload.vertice = inst.vertices.findAt( ((pload.x,
                                                pload.y,
                                                pload.z),) )

    if self.separate_load_steps == True:
        self.step2Name = 'Axial_Displ_Step'
        self.step1Name = 'PLoad_Step'
        # creating perturbation load step
        abaqus_functions.static_step(
            model                  = mod,
            name                   = self.step1Name,
            previous               = 'Initial',
            damping                = self.artificial_damping,
            initialInc             = 1.,
            maxNumInc              = self.maxNumInc,
            minInc                 = self.minInc,
            maxInc                 = self.maxInc,
            stabilizationMagnitude = self.stabilizationMagnitude,
            adaptiveDampingRatio   = self.adaptiveDampingRatio,
                                    )
        # creating axial load step
        abaqus_functions.static_step(
            model                  = mod,
            name                   = self.step2Name,
            previous               = 'PLoad_Step',
            damping                = self.artificial_damping,
            initialInc             = 0.1,
            maxNumInc              = self.maxNumInc,
            minInc                 = self.minInc,
            maxInc                 = self.maxInc,
            stabilizationMagnitude = self.stabilizationMagnitude,
            adaptiveDampingRatio   = self.adaptiveDampingRatio,
                                    )
    else:
        self.step1Name = 'PLoad_Axial_Displ_Step'
        self.step2Name = self.step1Name
        # creating perturbation load and axial displacement at same step
        abaqus_functions.static_step(
            model                  = mod,
            name                   = self.step2Name,
            previous               = 'Initial',
            damping                = self.artificial_damping,
            maxNumInc              = self.maxNumInc,
            initialInc             = 0.1,
            minInc                 = self.minInc,
            maxInc                 = self.maxInc,
            stabilizationMagnitude = self.stabilizationMagnitude,
            adaptiveDampingRatio   = self.adaptiveDampingRatio,
                                    )
    # applying boundary conditions
    localCsys = inst.datums[1]
    icount = 0
    for edge in self.inst_edges:
        icount += 1
        #bot boundary conditions
        abaqus_functions.displacementBC(\
                model          = mod,
                name           = 'BotBC_' + str(icount),
                createStepName = 'Initial',
                region         = Region( edges = edge.bot ),
                u1  = SET, u2  = SET, u3  = SET,
                ur1 = SET, ur2 = SET, ur3 = SET,
                localCsys      = localCsys )
        #top boundary conditions
        abaqus_functions.displacementBC(\
                model          = mod,
                name           = 'TopBC_' + str(icount),
                createStepName = 'Initial',
                region         = Region( edges = edge.top ),
                u1  = SET, u2  = SET, u3  = UNSET,
                ur1 = SET, ur2 = SET, ur3 = SET,
                localCsys      = localCsys )
        if self.axialInclude:
            abaqus_functions.displacementBC(\
                    model = mod,
                    name  = 'AxialDisplacement_' + str(icount),
                    createStepName = self.step2Name,
                    region = Region( edges = edge.top ),
                    u1  = UNSET, u2  = UNSET, u3  = -AXIALDISPL,
                    ur1 = UNSET, ur2 = UNSET, ur3 =       UNSET,
                    localCsys = localCsys )
    icount = 0
    if self.separate_load_steps == True:
        createStepName = 'PLoad_Step'
    else:
        createStepName = 'PLoad_Axial_Displ_Step'
    for pload in self.ploadconf.ploads:
        if pload.pltotal < TOL:
            print 'WARNING - pload less than TOL, will be ignored'
            continue
        icount += 1
        # creating the perturbation loads
        mod.ConcentratedForce( name='PLoad_' + str(icount),
                               createStepName = createStepName,
                               region= Region( vertices = pload.vertice ),
                               cf1 = pload.plx,
                               cf2 = pload.ply,
                               cf3 = pload.plz,
                             )
    # meshing the part
    part.generateMesh()
    # creating job
    job = mdb.Job( name =JOBNAME,
                   model = mod,
                   scratch = r'c:\Temp\abaqus\scratch',
                   memory = 4,
                   memoryUnits = GIGA_BYTES,
                   #numCpus = 6,
                 )
    self.job = job
    # defining output requests
    mod.fieldOutputRequests['F-Output-1'].setValues(\
                    variables = self.output_requests )

    if self.request_stress_output:
        print 'WARNING - stress outputs on - this will create a huge .odb file'
        mod.FieldOutputRequest(
            name                = 'stress_outputs',
            createStepName      = 'Axial_Displ_Step',
            variables           = ('S', 'HSNFTCRT', 'HSNFCCRT',
                                        'HSNMTCRT', 'HSNMCCRT'),
            layupNames          = ('InstanceCylinder.CompositePlate',),
            layupLocationMethod = ALL_LOCATIONS, 
            rebar               = EXCLUDE )
    # reading and mapping nodes
    self.read_and_map_nodes()

    # creating sets (will be used in the output recovery)
    
    for cross_section in self.cross_sections:
        nodes = []
        for node in cross_section.nodes:
            if node == None:
                continue
            nodes.append( node.obj )    
        mesh_node_array = mesh.MeshNodeArray( nodes=nodes )
        name = 'CROSS_SECTION_%03d' % cross_section.index
        part.Set( name=name, nodes=mesh_node_array )
    
    if self.imperfection_activate:
        # editing "Keyowords" to include the imperfection file
        text  = '** ----------------------------------------------------------'
        text += '\n**'
        text += '\n*IMPERFECTION, STEP=%d, FILE=%s' \
              % ( int(self.imperfection_step_num), self.imperfection_filename )
        list1 = self.imperfection_buckling_modes
        list2 = self.imperfection_scaling_factors
        for i in range(len( list1 )):
            text += '\n%d, %f' % ( int( list1[i] ), float( list2[i] ) )
        mod.keywordBlock.synchVersions( storeNodesAndElements = False )
        sieBlocks = mod.keywordBlock.sieBlocks
        pattern = '*Step'
        index = None
        for i in range(len( sieBlocks )):
            sieBlock = sieBlocks[ i ]
            if sieBlock.find( pattern ) > -1:
                index = i-1
                break
        if index == None:
            print 'WARNING - *IMPERFECTION entry was not created !'
            print '          *Step pattern not found !'
        else:
            mod.keywordBlock.replace( index, text )
        

